MASSACHUSETTS INSTITUTE OF TECHNOLOGY 
ARTIFICIAL INTELLIGENCE LABORATORY 



A.I. Memo No. 877 February. 1986 



THE OUTER SOLAR SYSTEM FOR 210 MILLION YEARS 



James H. Applegate 
Michael R. Douglas 

Yekta Gursel 
Gerald Jay Sussman 

Jack Wisdom 



Abstract 

We have used a special purpose computer to integrate the orbits of the outer five 
planets for 100 Myr into the future and 100 Myr into the past. The strongest features in the 
Fourier transforms of the orbital elements of the Jovian planets can be identified with the 
frequencies predicted by linear secular theory. Many of the weaker features in the Fourier 
spectra are identified as linear combinations of the basic frequencies. We note serious 
differences between our measurements and the predictions of Bretagnon (1974). The 
amplitude of the 3.796 Myr period libration of Pluto's longitude of perihelion is modulated 
with a period of 34 Myr. Very long periods, on the order of 137 million years, are also seen. 
The orbit of Pluto is stable for the duration of our integration; the maximum Lyapunov 
characteristic exponent is less than 10 " 6 ' 8 yr " l . 
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Abstract 

We have used a special purpose computer to integrate the orbits of the outer five 
planets for more than 100 Myr into the future and more than 100 Myr into the 
past. The strongest features in the Fourier transforms of the orbital elements of 
the Jovian planets can be identified with the frequencies predicted by linear secular 
theory. Many of the weaker features in the Fourier spectra are identified as linear 
combinations of the basic frequencies. We note serious differences between our 
measurements and the predictions of Bretagnon (1974). The amplitude of the 3.796 
Myr period libration of Pluto's longitude of perihelion is modulated with a period of 
34 Myr. Very long periods, on the order of 137 Myr, are also seen. The orbit of Pluto 
is stable for the duration of our integration; the maximum Lyapunov characteristic 
exponent is less than 10 yr. 

Introduction 

The determination of the stability of the solar system is one of the oldest 
problems in astronomy, but despite considerable attention there is no proof 
of the stability of the system. Systems with two degrees of freedom which are 
coupled sufficiently weakly have been shown to be stable for all time by the 
Kolmogorov-Arnold-Moser (KAM) theory; the stability of the solar system 
has been demonstrated under the conditions that the masses of the planets, 
their eccentricities and their inclinations are sufficiently small (Arnold 1961). 
The actual solar system, however, does not meet the stringent requirements 
of KAM theory. Certainly, the great age of the solar system demands a high 
level of stability, but weak instabilities may still be present. Wisdom (1982, 
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1983) has shown that planetary systems which appear to be stable may in fact 
be slightly unstable (have a small positive Lyapunov characteristic exponent), 
and that this slight instability can manifest itself in dramatic and relatively 
sudden changes in an orbit. The stability of the solar system should thus not 
be taken for granted. 

Pluto's orbit is unique among the planets. It is both highly eccentric 
(e « 0.25) and highly inclined (i « 16°). The orbits of Pluto and Neptune 
cross one another, a condition permitted only by the libration of a resonant 
argument associated with the 3/2 mean-motion commensurability. This res- 
onance assures that Pluto is at aphelion when Pluto and Neptune are in 
conjunction, thus preventing close encounters. While resonances may stabi- 
lize the short term evolution they also give rise to the most prominent zones 
of instability. The long-term stability of Pluto's orbit is thus of particular 
interest. 



Numerical integrations with the Orrery 

Historically, the long-term behavior of the solar system has been studied using 
perturbation theory. The interactions and the resulting motions are studied 
using series expansions which are expected to converge to the physically cor- 
rect answers. In the secular theory (the first step in all such investigations) 
only those terms representing averages that are independent of the longitudes 
of the planets are considered (Brouwer and van Woerkom 1950; Brouwer and 
Clemence 1961). This may be thought of as treating the planets as massive 
rotating rings instead of point masses orbiting the sun. Secular theory gives 
considerable insight into the long term evolution of the planets, but it is only 
approximate and rapidly becomes intractable as it is extended to higher orders 
(see Anolik et al. 1969; Bretagnon 1974, Duriez 1979). 

The development of computers made direct integrations of the equations 
of motion of the solar system possible. The first long-term integrations were 
made by Cohen and Hubbard (1965), who integrated the outer 5 planets for 
120,000 yr, and Cohen, Hubbard, and Oesterwinter (1973, hereafter CHO), 
who extended the duration of the integration to 1 Myr. The time span cov- 
ered by these integrations, though short compared to the age of the solar 
system, was long enough to reveal the previously unknown resonant stabiliza- 
tion mechanism for the orbit of Pluto described above (Cohen and Hubbard 
1965; Brouwer 1966), and to demonstrate (CHO) that the Brouwer and van 
Woerkom (1950) linear secular theory was a good approximation to many 
features of the motion of the outer five planets. Recently, the time interval 
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over which the integration has been performed was extended to 5 Myr by 
Kinoshita and Nakai (1984). 

Though direct, long-term numerical integrations have proven themselves 
to be powerful tools in celestial mechanics their usefulness is limited by the 
availability of computer resources. If performed on a DEC VAX 11/780, one 
of the 100 Myr integration runs that we report here would require roughly one 
year. The reason for the large computing requirement is simple. The solar 
system is 10 8 to 10 10 orbits old, depending on which orbit period is chosen for 
a reference, and the timestep in a direct integration must be a small fraction 
of the orbit period. The computing time may be reduced to some extent by 
integrating equations of motion with the highest frequency variations removed 
by averaging (e.g. Schubart, 1964; Williams and Benson, 1971; Wisdom, 1982). 
The averaging method has been quite useful for qualitative investigations. 
Unfortunately, the exact relationship between the averaged and unaveraged 
solutions is never clear. We integrate Newton's equations of motion for the 
solar system, modeling the bodies as point masses, as they stand. We have 
solved the computer power problem by designing and building our own special 
purpose computer, the Digital Orrery. 

The Digital Orrery is a specialized but programmable high-performance 
computer designed for the efficient numerical integration of the equations 
of motion for systems with a small number of bodies which move in roughly 
circular orbits, e.g. the solar system. Details of the design and construction of 
the Orrery are discussed in Applegate et al. (1985). In an jV-body system the 
interaction of a body with the N -I other bodies determines the force that will 
act on the body in each timestep. In a serial computer, the computation of the 
forces requires 0(N 2 ) time, while the integration step for N bodies requires 
O(N) time. The Orrery is a parallel computer with N concurrent processing 
units (planet computers) arranged in a ring. One processor is allocated for 
each body in the simplest programs. At each step every processor executes 
identical instructions on different data. This arrangement allows the force 
computation to be accomplished in O(N) time and the integration step to be 
performed in O(l) time. Additionally, the Orrery iV-body program is much 
more efficient than an TV-body program written for a conventional computer. 
There are two reasons for this: The Orrery iV-body program and the Orrery 
itself were designed together, with the result that there are essentially no 
wasted machine cycles; and the planet computers have a considerable amount 
of pipeline parallelism. Thus, though a planet computer does an isolated 64- 
bit floating point add or multiply in 1.25 microseconds (roughly as fast as a 
VAX 11-780 with a floating point accelerator) the Orrery can run a 10-body 
integration roughly 60 times faster than a VAX 1 1-780 with a floating point 
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accelerator. 

Our integrations 

We used the Orrery to perform several integrations of the outer planets ap- 
proximately 110 Myr forward in time and one integration approximately 107 
Myr backward in time. We used the CHO initial conditions and masses, with 
the exception that we ran several Plutos as zero mass test-particles with vari- 
ous initial conditions. Whenever we refer to "Pluto" we mean that particular 
test-particle with the CHO initial conditions for Pluto. The main analysis re- 
ported in this paper was performed on a composite of a forward 106776073.9 
yr run with a reverse 106776073.9 yr run for a total of 213552147.8 yr or 
roughly 214 Myr. 

We also did one run using the CHO mass of Pluto (Ms U n/(3-6 x 10 5 )). 
This mass is far too large, but this run allowed us to compare our results for 
the first 10 6 years with CHO. Over 110 Myr the results of this run do not 
differ greatly from the results of the test particle runs. 

We employ the 12th order Stormer predictor used by CHO with timesteps 
near 40 days in our outer planet integrations. Since the truncation error of 
our predictor is smaller than the roundoff error due to the finite word length 
of the machine no corrector is used. The Orrery has a 55-bit single precision 
mantissa (it has 64-bit single-precision floating point). The performance of the 
integrator was improved by using double precision in certain crucial operations 
(Applegate et al. 1985). Stormer's predictor is a linear multistep integrator 
for second order systems: 

x n+i = 2x n x n _i © h 2 (c a n + C!a n _i H h c m a n - m ) 

The circled operations are those where we used more than 64-bit precision. 

The integrations are started using a combination of a 4th order Runge- 
Kutta and the Stormer integrator. The Runge-Kutta is used to produce a 
starting set of accelerations and positions for the Stormer predictor. We start 
with a time step chosen to be less than a given value (we use 0.05 days) and 
such that the final time step is a power of two times the initial time step. 
The predictor is then used to advance the evolution until enough points are 
generated so that the Stormer predictor can be used again with double the 
step size. The step size is repeatedly doubled by this process until the desired 
step size is achieved. 

The numerical error in our calculation caused a linear increase in the to- 
tal energy of the system, as found by CHO, in contradiction to the prediction 
of Brouwer (1937) that the error in energy should grow as the square root 
of the number of steps. The linear energy error entails a quadratic longitude 
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error because of the energy dependence of the mean motion. The fractional 
energy error after 110 Myr is 2 x 10~ 8 . This corresponds to a fractional en- 
ergy error per year of 1.8 x 10~ 16 . The fractional energy error per year in 
CHO was 2.4 x 10~ 16 . By comparison, the relative mass lost by the Sun 
in 110 Myr due to electromagnetic radiation alone is about 7.4 x 10~ 6 (the 
mass lost due to solar wind is smaller by a factor of 4). Our energy error, 
which is almost entirely due to Jupiter, corresponds to an accumulated error 
in Jupiter's longitude of roughly 100 degrees in 100 Myr. This error will not 
affect the secular variations of the orbits; it is only significant for phenomena 
that depend upon the relative longitudes of the planets, such as mean-motion 
resonances. The solar system has a number of approximate commensurabil- 
ities but none are sufficiently exact for this integration error to change the 
low-order resonant structure of the system. The only exact commensurabil- 
ity is the Neptune-Pluto mean-motion resonance, but this resonance does not 
depend upon the longitude of Jupiter — the averaged integration of Williams 
and Benson (1971) shows it quite accurately. Thus the physical significance 
of much more accurate integrations is questionable. 

In the appendix we give the initial values of the system state (center-of- 
mass coordinates and velocities) we used for one typical 110 Myr integration 
and values of the system state that we obtained at 1 1 Myr intervals for the 
entire integration run. This should allow independent analysis of our inte- 
gration error. This run was done with a 40 day step size using CHO masses 
except that Pluto was given zero mass. 

Data and Interpretation 

Each run evolved the system for a total simulated time of about 110 Myr 
(about 10 9 timesteps with timesteps near 40 days). Such a run takes roughly 
one week of real time. The positions and velocities of the planets and the 
sun were read out once every 10 5 timesteps, giving us a sampling period near 
10,000 yr. The positions and velocities were converted to orbital elements, 
and the elements used to form the variables h, k, p and q defined as 

h = esin($"2 + u>); k = e cos(f2 + u>) 
p — sin(i/2) sin(n); q = sin(i/2) cos(17) 
where t is the inclination, e is the eccentricity, f2 is the longitude of the 
ascending node, and ui is the argument of the perihelion. Thus w — fi + u 
is the longitude of perihelion. The linear secular theory is particularly simple 
when written in terms of the variables h, k, p and q. The variables p and q 
specify the orientation of the orbital plane, and the variables h and k specify 
the eccentricity and the orientation of the orbit in the orbital plane. 
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Data in the time domain are difficult to digest, so most of the analy- 
sis was done in the frequency domain. High-resolution power spectra were 
obtained with the chirp z-transform algorithm (Rabiner, Schafer, and Rader 
1969). Spectra were examined using Hanning windows to suppress spectral 
leakage. We also made some low-resolution spectra with the ordinary fast 
Fourier transform (FFT) algorithm. 

Our data is sampled at rather long periods (about 10,000 years between 
samples) because of limited computer power available for analysis of the result- 
ing data (even at this rate a 110 Myr run yields about 5 Mbytes of ephemeris!). 
Consequently the spectra we obtain are not band-limited, so many of the lines 
are aliased versions of lines with periods shorter than 20,000 years. Fortu- 
nately, aliased lines can be identified because their frequencies change if the 
sampling rate is changed. Real lines are stable under such change of sampling 
rate. Thus we use the results of several runs with slightly different sampling 
periods to pick out the unaliased lines. These are then presented as tables of 
frequencies and amplitudes. Raw spectra of Jupiter's h and p are shown in 
Figure 1 and Figure 2 respectively. These figures show the full range from zero 
frequency to the Nyquist frequency, 1.25 x 10~ 7 cycles/day (about 1/20,000 
yr). The dominant line in Jupiter's h corresponds to a period of 305,505 yr 
and the large line in the center of the plot corresponds to a period of 45,906 
years. 

The frequencies and periods of the lines identified as the linear h — k and 
p- q features are listed in Table 1, along with the Brouwer and van Woerkom 
(1950) and Bretagnon (1974) predictions for the frequencies of these modes. 
In linear secular theory the interactions between the planets are averaged over 
mean longitudes, and only the lowest nonvanishing terms in eccentricity and 
inclination are kept. At this level of approximation the equations of motion 
are linear and h and k are decoupled from p and q, yielding two normal 
modes of the system for each planet in the system. The motion of each planet 
is composed of a combination of the normal modes. 

The conservation of angular momentum constrains the p - q equations 
so that one mode has zero frequency. Thus, for four planets there are 4 
independent normal modes for the h — k spectra, and 3 independent normal 
modes for the p — q spectra. The p — q frequencies are denoted g^, <77, and 
<7e; the h - k frequencies are denoted /s, /e, /V, and /g. g$ is the zero p - q 
frequency. The full iV-body system has 3N - 3 independent frequencies, one 
for each degree of freedom. Thus the four planet (five body) system has 12 
independent frequencies. Four of these are the mean motions, leaving the 8 
secular frequencies, one of which is zero. 

Nonlinearities in the secular theory produce linear combinations of the 
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Figure 1: Raw spectrum of Jupiter's h. The common logarithm of the amplitude is 
plotted against frequency (in cycles per day). The full range of frequencies, from 
to 1.25 x 10~ 7 day _1 is shown. The large line in the center of the plot is the / 6 line. 
The largest line on the left is the f$ line. 
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Figure 2: Raw spectrum of Jupiter's p. The common logarithm of the amplitude is 
plotted against frequency. The full range of frequencies, from to 1.25 x 10" day 
is shown. The large line in the center of the plot is the g 6 line. 
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fundamental frequencies. The terms of the secular disturbing function are all 
even in the inclinations and in the eccentricities. For each term the sum of the 
coefficients of the longitudes of perihelion and of the longitudes of ascending 
nodes must be zero so that the term is independent of the choice of reference 
direction. The sum of the coefficients of the longitudes of ascending nodes 
must be even so that the term is invariant under inversion of the reference 
plane. As a consequence, the sum of the coefficients of the longitudes of 
perihelia is also even. The amplitude of each term in the disturbing function is 
proportional to powers of the eccentricities, the inclinations, and the masses of 
the planets involved. In each term the difference of the power of the inclination 
of each of the planets and the absolute value of the coefficient of the longitude 
of the corresponding node is even. Similarly, the difference of the power 
of the eccentricity and the absolute value of the corresponding coefficient of 
the longitude of perihelion is even. Since for secular terms the sum of the 
coefficients of the longitudes of ascending node must be even, the sum of the 
powers of the inclinations must also be even; and likewise for the eccentricities, 
the sum of the powers of the eccentricities must be even. 

Short period terms, those which depend on the mean longitudes, may also 
contribute to the long period variations. In classical canonical perturbation 
theory these terms may be removed at any particular order of the masses 
by a von Zeipel transformation. The new canonical variables differ from the 
old variables by a quantity of one higher order in the masses than the term 
eliminated from the disturbing function but divided by the frequency of the 
argument of the term. The new secular Hamiltonian is obtained as an average 
of the Hamiltonian expressed in terms of these new variables over the mean 
longitudes. A term can only have a non-zero average if it is multiplied by 
a term with the same frequency. Consequently, the sum of the powers of 
the eccentricities and the sum of the powers of the inclinations of any new 
term in the secular Hamiltonian are both even. Thus secular terms which 
arise in this way satisfy the same constraints as the original terms in the 
secular Hamiltonian. Since such terms contain two factors of the frequency 
in the denominator they only have a significant contribution if the frequency 
is small, i.e. near resonant. For example, the Hill corrections included by 
Brouwer and van Woerkom are of this type. 

When written in terms of the variables h, k, p and q, the sum of the powers 
of h and k is even, and the sum of the powers of p and q is even. The equations 
of motion for h, k, p, and q contain derivatives of the disturbing function with 
respect to these variables. Thus the terms in the solutions for h, k, p, and q 
must be made up of odd numbers of fundamental frequencies. Indeed, terms 
in h and A; must have an odd number of / frequencies and an even number of 
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g frequencies; terms in p or q must have an even number of / frequencies and 
an odd number of g frequencies. Thus there should be no lines shared between 
the h - k and p - q solutions. Terms which produce combinations of three 
basic frequencies are two powers of small quantities smaller than the linear 
terms. Since the eccentricities, inclinations and masses are all small, linear 
combinations of small numbers of the basic frequencies should dominate. 

Table 2 lists major lines in the h — k and p-q spectra of the Jovian planets 
along with a proposed identification in terms of the fundamental frequencies 
identified with the frequencies predicted by the linear secular theory. Each 
entry in Table 2a represents a contribution to h and k: 

h = \_. a t sin(27rr,t + S t ) 

i 

k = y^ ati cos(2;rr,< + 5,); 

i 

each entry in Table 2b represents a contribution to p and q: 

P = ^2^s'm(2ns t t + e t ) 

i 

q = 2_. A cos(2tts,< + e,). 

i 

The amplitudes a, and 0i for each planet are represented by their common 
logarithms; the frequencies s, and r, are in inverse days. 

These tables were produced by the following process. We started with 
the rectangular coordinate data, sampled at about 10000 year intervals, for 
a 214 Myr composite run. The 20,000 samples were converted to orbital 
elements and then to h, k, p and q for each planet. (In Table 3 we provide 
a summary of the minimum, mean, and maximum values of the semimajor 
axes, the eccentricities, and the inclinations of the outer planets for a 214 
Myr run.) We used a chirp z-transform to make high-resolution spectra with 
200000 frequency-domain samples covering the full range of frequencies from 
zero to the Nyquist frequency. Possible lines were identified by finding all 
points in the spectrum higher than their nearest neighbors. The logarithms 
of the amplitudes at the peak were then fit with a parabola to determine the 
center frequency and amplitude of each potential line. The phases of each of 
the potential lines were determined by fitting a straight line to the phases of 
the same points. A typical h or k spectrum has about 1000 possible lines of 
logarithmic amplitude larger than -6.5. A typical p or q spectrum has about 
200 possible lines of logarithmic amplitude larger than -6.5. 

For each planet, the h and k line tables were compared. Potential lines 
were matched if the frequencies were within 5.0 x 10 _12 day _ and if their 
amplitudes matched to within 0.1. If such a match was found the line was 
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accepted with the average of the matching frequencies and amplitudes. The 
phases of both the h and k were preserved. Similar processing was done to 
merge the p and q spectra. Most lines passed this test. No potential lines 
were lost in h - k spectra with logarithmic amplitude greater than -5 or in 
the p - q spectra with logarithmic amplitude greater than -5.5. 

To remove aliases, the 214 Myr line tables were compared with line tables 
from two 110 Myr runs with different sampling rates. The 214 Myr composite 
run was sampled at intervals of 3999996 days. One 110 Myr comparison run 
(call it TP1) was sampled at intervals of 4000000 days, and the other 110 Myr 
comparison run (call it TP2) was sampled at intervals of 3873324 days. The 
TPl comparison line table was produced using the same chirp z-transform as 
was used in the 214 Myr run. The TP2 comparison line table was produced 
using an 8192 point Fourier transform. A line was considered real (unaliased) 
if it appeared in all three line tables, with a logarithmic amplitude tolerance 
of 0.1 and with a frequency tolerance of 5.0 x 10~ 12 day -1 in the TPl table 
and 1.0 x 10 _11 day _1 in the TP2 table. This reduced the number of distinct 
frequencies to 103. 

Even with high resolution spectra, the measurement of phase is numer- 
ically difficult. By contrast, the consistent measurement of frequencies and 
amplitudes was much easier; frequencies often matched predicted values to 
5 significant places. For each distinct unaliased h - k line we compared the 
h phase with the k phase (similarly for p - q) to determine the sign of the 
frequency (because these should lead or lag by jt/2). In a few cases the phases 
were inconsistent or ill defined. In these cases we were unable to give a good 
phase estimate or be sure of the sign of the frequency. No line of logarithmic 
amplitude larger than -5 had a poorly defined phase. In Table 2 we mark oc- 
currences of lines where there was substantial error in measurement of phase. 

We then associated lines of the same frequency that occurred in the 
corresponding spectra of more than one planet. An identification was made 
if the frequencies matched to within 5.0 x 10 _12 day _1 . We then checked 
the phases and the signs of the frequencies for consistency (phases had to 
be within 0.2 radians of each other, allowing a shift by n). This caused the 
disambiguation of a few lines, and allowed us to determine the phases of a 
few lines for which either the h or the k phase had serious measurement error. 
Consistent sets were then given averaged frequencies and phases. 

We hand selected the fundamental frequencies {/«,<?»} to be the largest 
lines in the vicinity of the Brouwer and van Woerkom predictions. We then 
generated all possible additive combinations of up to three fundamental fre- 
quencies. We attempted to identify each line with one of these combinations 
according to frequency and phase. Just as the frequencies of a combination 
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are formed by adding and subtracting frequencies of fundamental lines, the 
phases of a combination are formed by adding and subtracting the phases of 
fundamental lines. An identification was made if the frequencies matched to 
within 5.0 x 10 _12 day -1 and the phases matched to within 0.2 radians. These 
allowed further disambiguation of a few lines and confirmed the classification. 
In identification of lines there are 287 distinct combinations with 3 or fewer 
fundamental frequencies. Of the 103 lines we find, 23 could not be assigned 
to a combination of 3 frequencies. 

The amplitude coefficients may be either positive or negative. A negative 
amplitude may be viewed as a positive amplitude with the phase rotated by 
n. Since we are specifying the amplitudes with the logarithms of the absolute 
values, we also provide the sign. In Table 2 we indicate the sense of the match 
of the measured phase of a line with the phase implied by its identification 
label using the symbols "(+)" and "(-)". 

Even the comparison of only two runs eliminates most aliases. For the 
looser case (TP2) the spectra were effectively divided into about 240000 bins. 
There are 12000 bins in frequency and about 20 bins in amplitude (A/ ss 
1.0 x 10 _11 day _1 , f max - 1.25 x 10 -7 day~ x ). Thus the expected number of 
coincidences in a sample of 1000 randomly distributed frequencies is less than 
4. The use of two different runs of two different sampling rates essentially 
guarantees that no aliases survive. 

Comparison with the secular theory 

To date the most complete secular theory for the entire solar system is that of 
Bretagnon (1974). Except for the / frequency Bretagnon's predictions agree 
with our measurements more closely than do Brouwer and van Woerkom's 
predictions. Bretagnon's fa frequency appears to be in error by more than 
7 parts in 100. We assert with confidence that this difference is not a conse- 
quence of the fact that Bretagnon's theory is for eight planets and we only 
simulated the Jovian planets. To verify this fact we made a 3 Myr run of the 
Jovian planets with Venus, Earth and Mars. The fe frequency as determined 
in the 3 Myr run differed from the fa frequency as determined in our 1 10 Myr 
runs by only 1 part in 10 4 . The #6 frequency determinations differed by 3 
parts in 10 4 . 

Many of the lines we measure are identifiable with those predicted by 
Bretagnon. As expected, we observe no line that appears in both the h and 
p spectra of any planet. Each line is either an h line or a p line and it 
obeys the constraints on formation described above. The amplitudes of the 
fundamental frequencies are quite what is expected. We see none of the two- 
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frequency combinations, such as fs+95-97 (remember that# 5 = 0), predicted 
by Bretagnon. This is because our osculating elements are referenced to the 
invariable plane and Bretagnon uses ecliptic elements. We pick up those lines 
when we use ecliptic elements to represent the results of our integrations. 
Except for those lines, the spectra of h and p do not depend on the choice of 
the osculating element reference plane. 

However there are some major discrepancies between our measured am- 
plitudes and Bretagnon 's amplitudes for combinations of the fundamental 
frequencies. For example, the lines marked F and E in Table 2 cannot be con- 
structed as a combination of three or fewer fundamental frequencies. These 
lines are not included by Bretagnon, since he only included terms up to third 
order in inclinations and eccentricities and second order in the masses. How- 
ever, the amplitudes of these lines are considerably larger than the amplitudes 
of many of the terms included by Bretagnon. Line E in Saturn's h spectrum 
has an amplitude of 3.63 x 10 -4 (with logarithm -3.44). This is larger than 
all but 7 of Bretagnon's more than 200 corrections to the Lagrange solution 
for the Jovian planets. Even for identified lines there are serious differences. 
For example, line D is identified as fz - g$ + gi (in Bretagnon's notation 
V*5 _ 06 + #7)- Bretagnon predicts the amplitude for Jupiter to be 1.0 x 10~ 6 
and for Saturn to be 6.0 x 10~ 6 . We measure these to be 1.0 x 10 -4 and 
3.2 x 10~ 4 , respectively. The major discrepancies in h lines appear to be clus- 
tered around the 4.59 x 10 4 year fa line that we noted before to be discrepant 
in frequency, but this is not the whole story. Line T, in the p spectrum of 
Uranus is a low frequency line (period about 5.53 x 10 5 yrs.) which is predicted 
to be of amplitude 4.6 x 10 -5 but we measure it to be about 3.0 x 10 -4 . We 
note a large number of discrepancies between the measured amplitudes and 
the amplitudes predicted by Bretagnon. In some cases we obtain a different 
sign for the frequency of a line than predicted by Bretagnon. Bretagnon also 
predicts a number of very small lines that we do not detect. It is possible 
that these small differences are a result of measurement errors introduced by 
nearby larger lines in the spectrum. 

Some insight in the origin of these discrepancies can be derived by an 
examination of the fine structure of the spectrum near the f§ line. Figure 3 
shows a high-resolution view of this region of the raw spectrum of Jupiter's h. 
The large central line, labeled A, is the / 6 line. The / 6 line is surrounded by 
a cluster of smaller lines. The lines marked "*" are aliases of high-frequency 
undersampled components; they moved in runs with different sampling rates. 
We can see that the f$ line is bracketed by a symmetrical group of lines. 
The outermost members of the group, labeled B and C, are fs — fe — fi and 
/s + /e - fi- D and E, and F and G are also symmetrically placed around 
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Figure 3: Expanded view of Jupiter's h spectrum around the region of the / 6 line. The 
common logarithm of the amplitude is plotted against frequency. The line marked A 
is the f 6 line. Lines marked * are aliases. Note the symmetrical sextuplet of lines 
surrounding / 6 . The strong pair marked D and E are not taken into account in theories 
that only consider terms that are of third order in inclinations and eccentricities and 
second order in the masses. 
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/e- Evidently, the / 6 line is modulated with three dominant modulating 
frequencies. The line labeled D is identified as /s — g^ + g?. The symmetrical 
line, labeled E, has no identification. But because of the symmetry we can 
deduce that this line is actually -fs + 2fe + g$ - g-j. Similarly, since G is 
identified as fa - gs + gg we can deduce that F must be 2/ 6 - /V + <76 ~9s- Thus 
the major contribution to the amplitudes of the D and G lines must be fifth- 
order combinations rather than third-order combinations. For example, we 
may identify D as / 6 - (-/ 5 + fa + <j 6 - g 7 ) and E as fa + (-/ 5 + fa +9e -gi). 
This fifth-order labeling of E is confirmed by the agreement of the phase 
combination with the measured phase. The phase predicted by our labeling 
of the F line is not quite so accurately confirmed, but the F line is a very small 
line close to very large lines. The problem with Bretagnon's analysis is that 
higher-order terms that he did not include are more important than most of 
the terms he included. Duriez (1979) has constructed a theory (for only the 
outer planets) including seventh order terms, and indeed his prediction of the 
/e frequency is correct to 1 part in 300, but his f$ frequency is still in error 
by more than 1 part in 100. 

The Orbit of Pluto 

Pluto's orbit crosses that of Neptune and it is substantially inclined to the 
ecliptic. It is protected from a close encounter with Neptune by the libration 
of the resonant argument 3Api uto - 2ANeptune ~ E7piuto- This resonance has a 
period of about 19,857 years; the resonant argument varies about n between 
-1.5 and 1.5 radians (-86 to 86 degrees). Williams and Benson (1971) demon- 
strated that Pluto also takes part in a complex resonance of the argument of 
perihelion (w), with a period of about 3.796 Myr. Nacozy and Diehl (1978) 
developed a semi-analytical theory of this libration. Figure 4 shows the evo- 
lution of the argument of perihelion over the full 214 Myr. The 3.796 Myr 
libration about n/2 is amplitude modulated with a period near 34 Myr. 

In the power spectrum of the resonant argument there are numerous lines 
near that point in the spectrum where the 19,857 year libration period should 
appear. (Since the libration period is undersampled the group appears in an 
aliased position in the spectrum.) Figure 5 shows a expanded view of this 
region. What we see is a broad peak centered around a period of 19,857 years 
consisting of numerous sharper peaks separated by a frequency corresponding 
to a period of 3.69 Myr. This spectrum is characteristic of narrow-band 
frequency modulation. Evidently, the frequency of the mean-motion libration 
is modulated by the circulation of the longitude of the ascending node. 

The variation of the argument of perihelion is locked to the variations 



pi 







!!iiiij!!liliii!!!!iliiillll!iiii!l!!liiii!!J!!l!ii!li!i! 

^";:; 'i;'.-'.; 1 , :!,-;:,, :;':: ; !L:'|! i;:;iii!!!Pi!M : !! i ; !i ,: fi'ii;j;, , ^: ! :fJi i : = ;-;: 



; " if s r !/ U 



pl||iiii|pwi|lii||p!l| ( i"iijpj 



•4E10 4E10 

Figure 4: The argument of perihelion, w, of Pluto for 214 Myr. The abscissa is time 
measured in days. The 3.80 Myr libration is amplitude modulated with a 34 Myr 
period. 
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Figure 5: Expanded view of the power spectrum of the resonant argument 3Api uto - 
2A Nept un e - wpiuto of Pluto around the region where the 20000 year libration is 
aliased. In this case the aliased lines are reflected once through the Nyquist frequency, 
1.25 x 10~ 7 day _1 . The cluster of peaks is centered on the 19857 year line. The 
elements of the cluster are separated by the frequency corresponding to the 3.69 
Myr circulation of the ascending node. This pattern is characteristic of narrow-band 
frequency modulation. Note the unusually noisy appearance of this spectrum. 
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of the eccentricity and the inclination (Williams and Benson 1971). This 
coupling is illustrated in Figure 6 where the inclination is plotted versus the 
argument of perihelion. The 34 Myr variation appears as a slow expansion 
and contraction of the nearly circular trajectory on this i versus ui plot. The 
34 Myr period seen in the argument of perihelion is even more pronounced 
in the inclination. In Figure 7 we see that the 34 Myr modulation of the 
inclination has an additive as well as a multiplicative component. In addition 
there seems to be a much longer period present. This longer period is due 
to a near commensurability between <jg and the period of circulation of the 
ascending node. 

Because of the libration of the argument of perihelion the h - k and p - q 
degrees of freedom are more intimately coupled for Pluto than they are for 
the other planets. The h — k and p — q spectra of the Jovian planets have no 
lines in common; Pluto's h - k and p - q spectra have many lines in common. 
Now we expect the h — k and p - q spectra of Pluto to contain combinations 
of two new fundamental frequencies (which correspond to the two long period 
degrees of freedom of Pluto) as well as the fundamental frequencies associated 
with the motion of the Jovian planets. The two new fundamental frequencies 
must be chosen arbitrarily since there is no adequate long period theory of 
Pluto to serve as a guide; the semi-analytic theory of Nacozy and Diehl (1978) 
is over averaged and has only one degree of freedom. By analogy with the sec- 
ular perturbation theory of the other planets we choose the new fundamental 
frequencies from the h - k and p - q spectra (rather than, say, the spectra 
of the argument of perihelion). Figure 8 shows a greatly expanded view of 
the very low frequency portion of the h spectrum computed with the chirp 
z-transform algorithm. We choose the first fundamental frequency, p 1; to be 
the frequency of the circulation of the ascending node and the circulation of 
the longitude of perihelion. (These two frequencies are equal because of the 
libration of the argument of perihelion.) The largest peak in Figure 8, la- 
belled R, has the frequency pi . A suitable choice for the second fundamental 
frequency is less clear. We choose pi to be the frequency of the line marked 
W, which is the largest line in the long period portion of the h - k and p - q 
spectra which cannot be explained as a combination of the fundamental lines 
previously chosen. 

A summary of the main long-period lines and their relationships is given 
in Table 4. The rightmost peak in Pluto's h, labeled V, has a period near 
943,000 years. The largest peak, R, has the frequency p\, which corresponds 
to a period of 3.69 Myr. The second largest line, S, represents the long 
period that is visible in the time-domain h plot shown in Figure 9. The 
third largest line, T, is the gg line. Numerically, the frequency of line S 
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Figure 6: The inclination of Pluto to the invariable plane, on the ordinate, is plotted 
against the argument of perihelion, on the abscissa. Angles are measured in radians. 
The argument of perihelion librates about tt/2. The 34 Myr modulation produces the 
thickness of the band. 
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Figure 7: The inclination of Pluto for 214 Myr. Besides the 34 Myr modulation of the 
3.80 Myr oscillation, there is evidence for much longer period variations (or perhaps 
a secular drift!). 
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Figure 8: A highly expanded view of the low-frequency end of the power spectrum of 
Pluto's h. The main peak, R, is p x , the 3.69 Myr circulation of the ascending node 
and the longitude of perihelion. Line W is the p 2 line. The peak T is the g s line. 
Note the regular progression of peaks; the lines S, R, T, U, and V are all separated 
by the same frequency. Line S has a period near 137 Myr. 
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Figure 9: The h of Pluto for 214 Myr. The basic oscillation is the 3.69 Myr (pi) cir- 
culation of the ascending node and longitude of perihelion. The 27 Myr p 2 oscillation 
and the 137 Myr 2pi - g s oscillation are evident. The small cross indicates the value 
of h at the current epoch. 
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appears to be 2p\ - g&. Thus the frequency of the circulation of the node, 
Pi, is nearly commensurate with the fundamental frequency g$. Perhaps the 
commensurability accounts for the strength of the very long variations. This 
near commensurability could have significant consequences for the long-term 
stability of Pluto. 

The commensurability of g% with pi is an accidental property of Pluto's 
orbit. We have tried a variety of "Plutos" with nearby initial conditions. We 
have varied the initial inclination by ±5° without changing the other elements, 
and we have tried various values for the argument of perihelion between and 
?r/2, keeping the amplitude of the libration of the resonant argument con- 
stant by varying the initial semimajor axis. The observed commensurability 
between gg with pi does not exist for nearby initial conditions. 

The strongest components in the eccentricity should be simple differences 
of the frequencies of components of h — k. Inclination lines should be formed in 
a similar way from lines in p - q. In addition, frequencies of smaller lines may 
contain sums of frequencies of h - k or p - q lines and other even combinations 
of these lines. Indeed, the main lines in Pluto's e, i, and w are formed in this 
way. Table 5 shows the formation of several of the largest lines. 

An objective measure of the qualitative character of a solution of a set 
of differential equations is the Lyapunov characteristic exponent. A positive 
exponent indicates exponential divergence of neighboring trajectories; a zero 
exponent indicates quasiperiodic behavior. The maximum Lyapunov expo- 
nent can be estimated by computing the divergence of two nearby trajectories 
and plotting log7 = \og\\n(d(t) / d(to)) / (t - to)} versus log(< - to), where d(t) 
is the phase-space distance between the particles at time t and to is the ini- 
tial time. For chaotic trajectories log 7 will approach the logarithm of the 
Lyapunov exponent, while for quasiperiodic trajectories log 7 versus logt will 
follow a line with slope -1. In any finite computation we can only give an es- 
timate of this quantity which is, strictly speaking, only defined in the infinite 
time limit. Figure 10 shows the results of our computation of the Lyapunov 
exponent. Despite the odd long-period wanderings of Pluto's 1 and h, there 
is no objective evidence that Pluto, or the solar system as a whole, is chaotic. 
For a 110 Myr integration, the Lyapunov exponent of the solar system is es- 
timated to be less than 10 -6 8 /yr. We have experimented with some other 
test particles with initial conditions close to the initial conditions for Pluto. 
Some of these were clearly chaotic, yet they did not have a close encounter 
with another planet for a full 110 Myr integration. We will report on these 
in another paper. 
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Figure 10: The comparison of neighboring trajectories for Pluto shows no sign of a 
positive Lyapunov exponent. The absicssa is the common logarithm of the time since 
the start, in days, and the ordinate is the common logarithm of 7, measured in day - 
The maximum Lyapunov exponent is apparently less than 10" 
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Conclusions 

Our integrations support our prejudice that the solar system is quite stable 
over geologically significant time periods. 

Our integrations yield a model for the long-period variations of the ele- 
ments of the outer planets that is significantly more accurate than the best 
analytic solutions available to date. In particular, we find contributions to 
the variations that are not in Bretagnon's model which are larger than all 
but seven of his approximately 200 corrections to the Lagrange solution. Sur- 
prisingly, Brouwer and van Woerkom's estimate of the /g frequency is more 
accurate than Bretagnon's. By further integration experiments we have veri- 
fied that this discrepancy is not due to the effects of the inner planets. In fact, 
in the system of the Jovian planets alone, in the vicinity of the fa line there 
are large high-order contributions that are neglected in the analytic theory. 
This suggests that a very high order theory is necessary to get an accurate 
determination of the f& frequency by analytic means. 

We do not know the sensitivity of our results to small variations in the 
masses of the planets or the initial conditions. This is an ideal problem for 
investigation with a highly vectorized supercomputer where many solar system 
models could be run simultaneously. 

In any case, conclusions depending on the analytic theory must be re- 
considered in the light of our findings, especially those that depend on the 
exact frequency of the /q line. For example, the frequencies determine the 
positions of the secular resonances that in turn determine the positions of var- 
ious gaps in the asteroid belt and the inner edge of the asteroid belt. Because 
Williams and Faulkner (1981) use the Brouwer and van Woerkom frequencies 
to compute the secular resonance positions, their conclusions are not substan- 
tially modified — at zero inclination, the position of the /g secular resonance 
is shifted outward by only 0.02 AU from the position computed by Williams 
and Faulkner. 

The Milankovich hypothesis is that long-term climate variations on Earth 
are partially caused by variations in the insolation. The insolation is itself 
a function of the orbital elements of Earth and the orientation of its rota- 
tion axis. Berger (1976) has used Bretagnon's model in an evaluation of the 
Milankovich hypothesis. But a 7 percent difference in the value of the / 6 
frequency will cause the phase of that contribution to be inverted after only 
300000 years. On the other hand, the contribution of the fa line to the eccen- 
tricity of Earth's orbit is only about 1/6 of the largest contribution. In any 
case this work must be rechecked. 

We have observed extremely long period variations in the orbital ele- 
ments of Pluto. The eccentricity, inclination, and argument of perihelion of 
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Pluto are modulated with a 34 Myr period. Since the long period motion 
of Pluto has two degrees of freedom (corresponding to the two independent 
resonant arguments) the appearance of an additional fundamental frequency 
(p 2 at 27 Myr) is to be expected. The approximately 137 Myr variation we 
see in /ip] u to is more surprising. It indicates an accidental commensurability 
between the secular frequency, g&, and the frequency of the circulation of the 
node, pi. Perhaps this is important. While chaotic behavior of Pluto could 
have occurred in the two degree of freedom long-period problem, on the time 
scale of our integrations, we see no evidence of such instability. The maxi- 
mum Lyapunov characteristic exponent appears to be less than 10~ 6 yr -1 . 
However, with the 20000 year mean-motion resonance and the 3.796 Myr w 
libration resonance, this commensurability between g$ and p\ may provide 
just what is needed to make the orbit of Pluto chaotic over much longer time 
scales. Perhaps we will observe instability in our billion year integrations. 
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Tables 

Table la: Fundamental frequencies as measured, compared with the theoret- 
ical predictions of Brouwer and van Woerkom and of Bretagnon. The mea- 
sured frequencies are averages over the occurrences of the lines in the chirp 
z-transform spectra of each of the planets. Numbers in parentheses indicate 
the formal error (standard deviation) in the last digit shown. 





Measured 




Br-vW 


Bretagnon 




Period 


Frequency 




Frequency 


Frequency 




yr 


cycles/day 




cycles/day 


cycles/day 


18 


1 . 924950e6 


1.4222969e-9 


(41) 


1.33790e-9 


1.410887e-9 


17 


4.193511e5 


6.5287787e-9 


(51) 


5.74464e-9 


6.475315e-9 


16 


4.590574e4 


5.9640698e-8 


(11) 


5.86738e-8 


5.538393e-8 


f5 


3.055055e5 


8.9617245e-9 


(1) 


9.07527e-9 


8.887885e-9 


g8 


1.872103e6 


-1.4624465e-9 


(16) 


-1.43129e-9 


-1.460675e-9 


§7 


4.327497e5 


-6.3266380e-9 


(1) 


-6.13199e-9 


-6.337273e-9 


g6 


4.921700e4 


-5.5628146e-8 


(1) 


-5.43631e-8 


-5.549021e-8 


g5 

















Table lb: Phases of fundamental lines as measured, compared with the the- 
oretical predictions of Brouwer and van Woerkom and of Bretagnon. The 
measured phases are averages over the occurrences of the lines in the chirp 
z-transform spectra of each of the planets. Numbers in parentheses indicate 
the formal error (standard deviation) in the last digit shown. 





Measured 


Br-vW 


Bretagnon 




Radians 


Degrees 


Degrees 


Degrees 


f8 


1.2870 (16) 


73.742 


69.431 


72.090 


17 


2.0769 (5) 


118.997 


131.944 


114.775 


16 


2.2193 (12) 


127.155 


131.686 


127.715 


15 


0.5242 (2) 


30.032 


31.174 


28 . 503 


g8 


3.5521 (12) 


203.518 


202 . 293 


201.288 


g7 


5.5886 (1) 


320.202 


315.063 


316.293 


g6 


2.2228 (4) 


127.359 


127.366 


125.643 


g5 
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Table 2a: Dominant h - k spectral lines for the Jovian planets in order of increasing frequency. 
Amplitudes are represented by their common logarithms. Frequencies and phases are averages of 
measured components in chirp z-transforms. Phases are computed for J.D. 2430000.5. Components 
marked ( + ) have the given phase; components marked (-) have phase shifted by n from the given 
phase. 



Identity 


Period 
years 


Frequency 
l/days 


Phase 
radians 


Jupiter 


Saturn 


Uranus 


Neptune 


-16+17+18 


2.709005e+6 


-1.010648e-9 


2.8378 


--- 


--- 


--- 


-4.26 (+) 


+18 


1.924960e+6 


1.422297e-9 


1 . 2870 


-4.23 (+) 


-4.23 (+) 


-2.77 (+) 


-2.04 (+) 


+17+g7-g8 


1.64468U+6 


1.664670e-9 


4.1099 


... 


... 


--- 


-3.81 (-) 


+18+g7-g8 


7.964332e+6 


-3.441962e-9 


3.3854 


... 


... 


... 


-6.64 (-) 


-17+18+18 


7.431357e+6 


-3.684187e-9 


0.4977 


... 


... 


... 


-6.43 (+) 


+15-17+18 


7.101583e+6 


3.866268e-9 


6.0136 


-6.60 (-) 


... 


-4.30 (-) 


-4.31 (-) 


-16+17+17 


6.684161e+6 


4.096027e-9 


3.6063 


-6.00 (+) 


... 


... 


... 




6.682436e+6 


4.097086e-9 


1.7286 


... 


... 


... 


-4.84 (+) 




6.681535e+6 


4.097636e-9 


3.4629 


... 


... 


-3.19 (+) 


... 


+16+g6-g8 


5.000601e+6 


5.475163e-9 


0.9037 


... 


... 


-4.76 (+) 


... 


-16+18+18 


4.476731e+6 


-6.117103e-9 


2.0439 


... 


... 


... 


-6.68 (-) 


+18-g7+g8 


4.364766e+6 


6.287022e-9 


6.6236 


-6.36 (+) 


... 


-4.12 (+) 


... 


+17 


4.193611e+6 


6.628779e-9 


2.0769 


-2.73 (+) 


-2.81 (+) 


-1.64 (+) 


-2.43 (-) 




3.609930e+5 


7.684221e-9 


2 . 1969 


-6.36 (+) 


... 


... 


--- 




3.462663e+6 


?7.906777e-9 


?? 


... 


--- 


-6.47 


--- 


+16 


3.066049e+6 


8.961726e-9 


0.6242 


-1.36 (+) 


-1.48 (+) 


-1.42 (-) 


-2.71 (+) 


+16+g6-g7 


2.647994e+6 


1.033934e-8 


6.1811 


-6.07 (+) 


-4.96 (+) 


-4.31 (-) 


-6.55 (+) 


+16+16-17 


2.402761e+6 


1.139466e-8 


6.2796 


-4.16 (-) 


-4.34 (-) 


-3.34 (+) 


-4.20 (-) 


+17+17-18 


2.363096e+6 


1.163610e-8 


2.9233 


... 


... 


-4.98 (-) 


-5.13 (+) 


+16-g7+g8 


1.980214e+5 


1.382603e-8 


4.6633 


--- 


--- 


... 


-5.26 (+) 




1.980122e+6 


1.382668e-8 


1.2765 


... 


... 


-5.10 (+) 


... 


+16+17-18 


1.946128e+6 


1.406820e-8 


1.3116 


... 


... 


-4.89 (+) 


-4.86 (-) 


-17+g7+g8 


1.912215e+6 


-1.431769e-8 


0.7791 


--- 


--- 


-6.36 (+) 


-6.64 (-) 


+15+15-18 


1.669193e+5 


1.660110e-8 


6.0467 


-6.08 (-) 


--- 


... 


-6.15 (+) 


-16+g7+g8 


1.634486e+6 


-1.676063e-8 


2.2897 


... 


--- 


-6.36 (-) 


-6.67 (+) 


-17+g7+g7 


1.427337e+6 


-1.918163e-8 


2.7827 


... 


... 


-6.39 (-) 


-5.80 (+) 




1.266758e+6 


?2.161305e-8 


?? 


... 


--- 


-6.07 


--- 


+16+g6-g7 


6.787014e+4 


-4.033965e-8 


3.4610 


... 


... 


-6.22 (+) 


... 


+16+16-16 


6.662876e+4 


-4.171726e-8 


6.1131 


-4.86 (+) 


-3.86 (-) 


-4.37 (+) 


... 


+16-16+17 


6.201225e+4 


-4.415016e-8 


0.3790 


-6.48 (+) 


-4.74 (-) 


-4.42 (-) 


... 


-16+16+18 


5.254864e+4 


6.210127e-8 


2.9811 


-6.46 (+) 


-4.98 (-) 


... 


--- 


+16-17+18 


5.020428e+4 


6.463421e-8 


1.4081 


-6.08 (+) 


... 


... 


... 


+17-g6+g7 


4.903883e+4 


6.683026e-8 


6.4446 


-6.32 (+) 


-4.84 (-) 


... 


— 




4.876614e+4 


6.616397e-8 


2 . 5887 


-6.48 (+) 


... 


... 


... 


B -16+16+17 


4.786814e+4 


6.720763e-8 


3.7831 


-3.70 (+) 


-3.23 (-) 


... 


... 


D +16-g6+g7 


4.699110e+4 


6.826317e-8 


3.8835 


-3.98 (-) 


-3.60 (+) 


-4.74 (-) 


... 


F 


4.673179e+4 


6.868647e-8 


1.3468 


-6.07 (-) 


-4.69 (+) 


... 


... 


A +16 


4.690676e+4 


6.964070e-8 


2.2193 


-1.80 (-) 


-1.32 (+) 


-2.81 (-) 


-3.99 (-) 


G +17-g6+g8 


4.510871e+4 


6.069461e-8 


3.3787 


-6.13 (+) 


... 


--- 


--- 


E 


4.486943e+4 


6.101818e-8 


0.6474 


-3.93 (+) 


-3.44 (-) 


... 


... 
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C +16+16-17 
+16-g6+g8 



+16+17-18 
+16+16-18 



+16-g6+g7 

-15+16+16 

+16+16-17 
+16+16-18 



4 


410649e+4 


6 


207366e-8 





6658 


-3.71 


(-) 


4 


337023e+4 


6 


312742e-8 


1 


8581 


-6.40 


(-) 


4 


314901e+4 


6 


345107e-8 


2 


1368 


-6.47 


( + ) 


4 


244300e+4 


6 


460663e-8 


2 


1316 


-5.77 


( + ) 


4 


228564e+4 


6 


4746B9e-8 


2 


9744 


--- 




4 


075388e+4 


6 


718013e-8 


1 


4533 


-6.48 


(-) 


3 


860693e+4 


-7 


091788e-8 


3 


5400 


--- 




2 


570631e+4 


1 


065091e-7 


4 


0056 


--- 




2 


537709e+4 


1 


078867e-7 


E 


4694 


-4.82 


( + ) 


2 


513122e+4 


1 


089422e-7 


E 


6879 


-6.14 


(-) 


2 


48206U+4 


?1 


103066e-7 




?? 


--- 




2 


481743e+4 


1 


103197e-7 


3 


9121 


-3.24 


(-) 


2 


4B1136e+4 


1 


116972e-7 


2 


2469 


-6.07 


( + ) 


2 


428192e+4 


1 


127B27e-7 


2 


3600 


-4.30 


(-) 


2 


322986e+4 


1 


178591e-7 


3 


1472 


-5.71 


(-) 



3 


21 


( + ) 


4 


90 


( + ) 


4 


98 


(-) 


5 


28 


(-) 


5 


66 


( + ) 


4 


99 


( + ) 


5 


84 


( + ) 


4 


29 


(-) 


4 


60 


( + ) 


2 


72 


( + ) 


4 


64 


(-) 


3 


77 


(+) 


- 


-- 





-6.70 (+) 



-4.86 (+) 
-4.38 (-) 



-5.61 (-) 



Phase measurement errors: 

? in sign of frequency: sign of frequency undetermined - bad phase 
? after phase: best phase estimate assuming positive frequency 
?? in place of phase: phase errors too large to determine phase 
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Identity 


Period 


Frequency 




years 


1/days 


+16-f7+g8 


2.821078e+6 


9.704982-10 


+f7-18+g7 


2.244632e+6 


-1.219787e-9 


+ g8 


1.872103e+6 


-1.462447e-9 


+16-17+g6 


1.088082e+6 


-2.616218e-9 


+17-i8+g8 


7.613263e+6 


3.644029e-9 


+I6-i7+g7 


7.031687e+5 


-3.893590e-9 




7.029241e+5 


-3.894946e-9 




6.619870e+B 


-4.136808e-9 


-16+i6+g6 


6.631944e+6 


-4.949166e-9 




4.665939e+6 


-6.009411e-9 


+f6-f8+g8 


4.506196e+6 


?6.077098e-9 


+ g7 


4.327497e+6 


-6.326638e-9 


-17+i8+g8 


4.168010e+5 


-6.668724e-9 




3.708698e+6 


-7.382243e-9 


-15+f7+g7 


3.125650e+6 


-8.759681e-9 


-f5+18+g8 


3.041431e+6 


-9.001860e-9 




2.446399e+5 


-1.119135e-8 




2.446353e+6 


-1.119166e-8 


-17+i8+g7 


2.394676e+5 


-1.143308e-8 


+15+i8-g8 


2.311113e+5 


1.184646e-8 


-16+18+g7 


1.974497e+5 


-1.386607e-8 


+17+f8-g7 


1.917571e+5 


1.427770e-8 


+17+i7-g8 


1.885576e+6 


1.461997e-8 


+f6+18-g7 


1.638386e+5 


1.671067e-8 


+*6+*7-g8 


1.614973e+5 


1.696292e-8 


+17+i7-g7 


1.412412e+5 


1.938423e-8 


-16-i5+g8 


1.412312e+5 


-1.938569e-8 


+16+17-g7 


1.264908e+6 


2.181716e-8 




1 . 196968e+6 


2.287322e-8 


+i6+*6-g7 


1.129007e+5 


2.425008e-8 




1. 026066* +6 


2.668298e-8 


-16+i6+g7 


6.172963e+4 


4.435229e-8 




5.987014e+4 


4.672982e-8 


+16-17+g7 


6.851960e+4 


4.678520e-8 


-15+i6+g8 


6.662908e+4 


4.921618e-8 


+i5-i6+g8 


5.260838e+4 


-5.214121e-8 


+16-f7+g6 


6.146803e+4 


-5.319617e-8 




5.046667e+4 


-6.425067e-8 


-16+i7+g8 


5.016749e+4 


-6.467420e-8 


+ g 6 


4.921701e+4 


-6.662816e-8 



Jupiter 



Saturn 



Uranus 



Neptune 



Table 2b: Dominant p - q spectral lines for the Jovian planets in order of increasing frequency. 
Amplitudes are represented by their common logarithms. Frequencies and phases are averages of 
measured components in chirp z-transforms. Phases are computed for J.D. 2430000.5. Components 
marked ( + ) have the given phase; components marked (-) have phase shifted by n from the given 
phase. 

Phase 

radians 

2.0022 

0.1404 

3.6521 

2.3641 

4.3627 

3.8723 

1.6392 

1.1958 

3.9167 

2.9019 

2.7833? 

6.6886 

2.7633 

2.3494 

0.8677 

4.3144 

4.8619 

2.0267 

4.7639 

4.6730 

. 0676 

4.0761 

0.5989 

2.6067 

5.3656 

4.9353 

2 . 6706 

3.3018 

3.4271 

1.7432 

3.3569 

1.0016 

6.6122 

5.7446 

5.3260 

1.8790 

0.6683 

3.6929 

3.2468 

2.2228 



-6.84 (-) 


-5.82 (-) 


-6.41 (-) 


-4.75 (+) 


-6.26 (+) 


-5.27 (+) 


-6.61 (-) 


-4.24 (-) 


-3.24 (-) 


-3.25 (-) 


-3.26 (+) 


-2.24 (+) 


— 


... 


-4.58 (-) 


-4.78 (+) 


... 


... 


... 


-5.82 (+) 


-4.96 (-) 


-5.00 (-) 


-3.74 (+) 


... 


... 


... 


... 


-4.70 (+) 


— 


... 


-6.86 (♦) 


— 


-4.76 (-) 


-6.23 (+) 


-3.53 (+) 


-4.25 (-) 


— 


... 


-6.16 (+) 


... 


... 


... 


... 


-6.01 (-) 


-3.32 (-) 


-3.41 (-) 


-2.06 (+) 


-2.97 (-) 


-5.61 (-) 


-6.61 (-) 


-4.23 (+) 


-6.07 (-) 


... 


... 


-6.48 (+) 


-6.03 (-) 


-5.00 (+) 


-6.15 (+) 


-3.71 (-) 


-4.66 (+) 


— 


... 


-6.10 (-) 


-6.73 (+) 


— 


... 


... 


-6.02 (+) 


— 


... 


-6.60 (+) 


... 


... 


... 


-6.58 (-) 


-6.76 (+) 


... 


... 


-6.16 (-) 


... 


— 


... 


-5.93 (+) 


-6.78 (-) 


— 


... 


-6.15 (-) 


-6.16 (+) 


— 


... 


-5.76 (-) 


-6.80 (+) 


... 


... 


... 


-6.11 (-) 


... 


... 


-6.45 (+) 


-6.46 (-) 


— 


... 


-6.52 (+) 


... 


— 


... 




-6.66 (+) 


— 


... 


-4.91 (-) 


-5.28 (+) 


... 


... 


-6.16 (+) 


... 


-6.08 (-) 


-5.94 (-) 


-4.71 (+) 


-5.61 (-) 


... 


... 


-6.12 (+) 


... 


... 


-5.89 (+) 


-6.12 (-) 


... 


... 


-6.09 (+) 


... 


... 





... 


-5.62 (+) 


... 


... 


... 


-6.13 (-) 


... 


-6.81 (-) 


-5.43 (+) 


-6.16 (+) 


... 


-6.42 (+) 


-4.98 (-) 


-6.33 (+) 


... 


-6.86 (+) 


-6.47 (-) 


... 


... 


-6.00 (+) 


-6.61 (-) 


... 


... 


-2.50 (-) 


-2.10 (+) 


-3.45 (-) 


-4.42 (-) 
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-16+17+g6 
-16+i7+g7 
+*6+*7-g6 
+i5+f6-g6 
+i6+f6-g7 
+15-i6+g6 
-f6+f7+g6 
+15+16-g6 



4.802774e+4 
4.716469e+4 
4.606187e+4 
3.849691e+4 
3.722364e+4 
3.663929e+4 
2.57B416e+4 
2.517795e+4 
2.203846e+4 



6 


700561e-8 


3 


8941 


6 


806104e-8 


3 


7830 


5 


9438B6e-8 


6 


4469 


7 


111871e-8 





3776 


7 


355159e-8 


6 


1076 


7 


492896e-8 


3 


4386 


1 


06307 le -7 





5274 


1 


087400e-7 


2 


0771 


1 


242306e-7 





5211 



-4.63 (+) 

-5.34 (-) 

-6.62 (-) 

-5.21 (-) 

-5.14 (+) 

-5.07 (+) 



4 


22 


(-) 


4 


98 


< + ) 


5 


18 


( + ) 


5 


99 


(-0 


4 


80 


( + > 


5 


81 


(-) 


5 


10 


(-) 


4 


70 


(-) 



4 


98 


( + ) 


6 


46 


( + ) 


6 


49 


(-) 


6 


81 


( + ) 


6 


64 


(-) 


6 


17 


( + ) 


6 


43 


(-) 


5 


88 


( + ) 


5 


33 


(-) 



Phase measurement errors: 

? in sign of frequency: sign of frequency undetermined - bad phase 
? after phase: best phase estimate assuming positive frequency 
?? in place of phase: phase errors too large to determine phase 
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Table 2c: Given the frequencies and phases measured for the fundamental lines we 
compute the frequencies and phases of the combination lines that occur in our data. 
Here we have the predicted (not measured) frequencies and phases, and the differences 
between the measured and predicted values. 





h - 


k lines 






Identity 


Frequency 


Error 


Phase 


Error 


-16+17+18 


-1.010649d-09 


9.000d-16 


2 . 8398 


-0.0020 


+18 


1.422297d-09 


1.000d-16 


1 . 2870 


-0.0000 


+17+g7-g8 


1.664687d-09 


8.280d-14 


4.1134 


-0.0036 


+18+g7-g8 


-3.441896d-09 


-6.740d-14 


3.3236 


0.0618 


-17+18+18 


-3.684186d-09 


-2.100d-15 


0.4972 


0.0006 


+15-17+18 


3.855243d-09 


2.630d-14 


6.0175 


-0.0039 


-16+17+17 


4.095833d-09 


1.941d-13 


3.6296 


-0.1243 


+16+g6-g8 


6.474998d-09 


1.546d-13 


0.8901 


0.0136 


-16+18+18 


-6.117131d-09 


2.770d-14 


2.0499 


-0.0060 


+18-g7+g8 


6.286488d-09 


6.336d-13 


6.6337 


-0.0101 


+17 


6.528779d-09 


3.000d-16 


2.0769 


0.0000 


+16 


8.961726d-09 


6.000d-16 


0.5242 


0.0000 


+16+g6-g7 


1.033919d-08 


1.600d-13 


6 . 1367 


0.0444 


+16+16-17 


1 . 139467d-08 


-2.030d-14 


6.2646 


0.0250 


+17+17-18 


1.163626d-08 


-1.606d-13 


2.8667 


0.0566 


+16-g7+g8 


1.382692d-08 


1.140d-13 


4.7708 


-0.1076 


+16+17-18 


1.406821d-08 


-6.300d-16 


1.3140 


-0.0024 


-17+g7+g8 


-1.431786d-08 


1.732d-13 


0.7806 


-0.0016 


+16+15-18 


1.660115d-08 


-5.210d-14 


6.0445 


0.0022 


-16+g7+g8 


-1.675081d-08 


2.790d-13 


2 . 3333 


-0.0436 


-17+g7+g7 


-1.918205d-08 


6.247d-13 


2.8171 


-0.0344 


+15+g6-g7 


-4.033978d-08 


2.335d-13 


3.4416 


0.0194 


+15+16-16 


-4.171726d-08 


-1.000d-16 


5.1122 


0.0009 


+15-16+17 


-4.416019d-08 


3.480d-14 


0.3818 


-0.0028 


-16+16+18 


5.210127d-08 


-4.000d-16 


2.9822 


-0.0011 


+16-17+18 


6.463422d-08 


-6.200d-16 


1.4294 


-0.0213 


+17-g6+g7 


6.683029d-08 


-2.670d-14 


6.4426 


0.0020 


-16+16+17 


6.720776d-08 


-1.222d-13 


3.7720 


0.0111 


+I5-g6+g7 


6.826323d-08 


-6.250d-14 


3.8899 


-0.0064 


+16 


6.964070d-08 


2.000d-16 


2.2193 


0.0000 


+17-g6+g8 


6.069448d-08 


3.180d-14 


3.4061 


-0.0274 


+16+16-17 


6.207364d-08 


6.200d-16 


0.6666 


-0.0007 


+16-g6+g8 


6.312742d-08 


-4.000d-16 


1.8634 


0.0047 


+16+17-18 


6.474718d-08 


-5.898d-13 


3.0091 


-0.0347 


+16+16-18 


6.718013d-08 


4.400d-16 


1.4564 


-0.0031 


+16-g6+g7 


1.089422d-07 


-6.000d-16 


5 . 6850 


0.0029 


-15+16+16 


1.103197d-07 


2.860d-14 


3.9144 


-0.0023 


+16+16-17 


1.127626d-07 


8.270d-14 


2.3617 


-0.0017 


+16+16-18 


1.178591d-07 


9.000d-16 


3.1615 


-0.0043 
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P - 


q lines 






Identity 


Frequency 


Error 


Phase 


Error 


+f5-f7+g8 


9. 704993d- 10 


-l.lOOd-lB 


1.9993 


0.0029 


+17-18+g7 


-1.220166d-09 


3.692d-13 


0.09B2 


0.0462 


+ g« 


-1.462447d-09 


-6.000d-16 


3.6621 


0.0000 


+*6-f7+g6 


-2.616227d-09 


8.700d-16 


2.3662 


-0.0011 


+17-f8+g8 


3.644036d-09 


-6.300d-16 


4.3419 


0.0108 


+iB-17+g7 


-3.893692d-09 


1.022d-13 


4.0368 


-0.163B 


-16+f6+g6 


-4.949173d-09 


7.B00d-lB 


3.9180 


-0.0013 


+f6-18+g8 


6.076981d-09 


1.169d-13 


2.7892 


-0.0069 


+ 87 


-6.326638d-09 


-1.034d-26 


6 . 6886 


0.0000 


-17+i8+g8 


-6.668928d-09 


2. 043d- 13 


2.7622 


0.0011 


-16+i7+g7 


-8.769684d-09 


2.800d-lB 


0.8681 


-0.0004 


-iB+18+g8 


-9.001874d-09 


2.410d-14 


4.3149 


-0.0006 


-17+i8+g7 


-1.143312d-08 


3.980d-14 


4.7987 


-0.0448 


+16+f8-g8 


1 . 184647d-08 


-7.900d-lB 


4.6423 


0.0307 


-*6+f8+g7 


-1.386607d-08 


-4.400d-16 


0.0683 


-0.0008 


+17+i8-g7 


1.427771d-08 


-1.360d-14 


4.0B86 


0.0166 


+i7+17-g8 


1.4B2000d-08 


-3.390d-14 


0.6017 


-0.0028 


+16+18-g7 


1.671066d-08 


1.060d-14 


2.6068 


0.0009 


+*5+f7-g8 


1.69B296d-08 


-2.970d-14 


6.3322 


0.0333 


+17+f7-g7 


1.938420d-08 


3.460d-14 


4.8484 


0.0869 


-fE-f6+g8 


-1.938690d-08 


3.0BBd-13 


2 . 6037 


0.0669 


+15+i7-g7 


2.181714d-08 


8.800d-16 


3.2967 


0.0061 


+16+f6-g7 


2.426009d-08 


-7.000d-lB 


1.7429 


0.0003 


-16+f6+g7 


4.43E234d-08 


-4.BB0d-14 


1.0006 


0.0011 


+16-17+g7 


4.678B28d-08 


-8.130d-14 


6.7310 


0.0136 


-f6+*6+g8 


4.921663d-08 


-3.470d-13 


5.2472 


0.0788 


+16-16+g8 


-6.214142d-08 


2.100d-13 


1.8669 


0.0221 


+16-i7+g6 


-6.319520d-08 


3.020d-14 


0.6701 


-0.0018 


-i6+17+g8 


-B.4B7437d-08 


1.6B8d-13 


3.4097 


-0.1639 


+ g 6 


-6.662816d-08 


-4.000d-16 


2.2228 


-0.0000 


+16-16+g7 


-6.700661d-08 


l.B00d-16 


3.8934 


0.0007 


-iB+i7+g6 


-B.806109d-08 


6.180d-14 


3.7766 


0.0074 


-f6+17+g7 


-B.9438B6d-08 


-2.700d-16 


6.4462 


-0.0003 


+16+f7-g6 


7.111866d-08 


6.080d-14 


0.3782 


-0.0007 


+lB+iB-g6 


7.366160d-08 


-6.000d-lB 


6.1087 


-0.0011 


+16+f6-g7 


7.492906d-08 


-1.006d-13 


3.4380 


0.0006 


+16-i6+g6 


-1.063071d-07 


1.960d-14 


. B277 


-0.0003 


-16+i7+g6 


-1.087401d-07 


6.B30d-14 


2.0804 


-0.0033 


+!B+i6-g6 


1 . 242306d-07 


3.160d-14 


. 6206 


0.0006 
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Table 3: Variations of the semimajor axes, a, the eccentricities, e, and the inclinations, 
t, of the outer planets over 214 Myr. a is in AU. i is in radians, i is measured relative 
to the invariable plane of CHO. 

mean 



minimum 

Jupiter 

a: 5.20135E+000 
e: 2.52206E-002 
i: 2.83801E-003 

Saturn 

a: 9.51297E+000 
e: 8.01269E-003 
i: 1.24408E-002 

Uranus 

a: 1.91016E+001 
e: 7.01759E-004 
i: 1.40198E-002 

Neptune 

a: 2.99143E+001 
e: 7.58230E-005 
i: 7.89724E-003 

Pluto 

a: 3.89838E+001 
e: 2.07307E-001 
i: 2.55640E-001 



5 . 20257E+000 
4.56015E-002 
6.45747E-003 

9.55488E+000 
5.36152E-002 
1 . 57810E-002 

1.92183E+001 
4 . 39952E-002 
1.78314E-002 

3.01098E+001 
1.01261E-002 
1 . 17275E-002 

3 . 94559E+001 
2.44315E-001 
2.78190E-001 



maximum 

5 . 20487E+000 
6.18055E-002 
9.76233E-003 

9 . 59206E+000 
8.90992E-002 
1.90587E-002 

1.93320E+001 
7.61826E-002 
2.15983E-002 

3.03168E+001 
2.28723E-002 
1.52102E-002 

3.99831E+001 
2.80581E-001 
2 . 97766E-001 
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Table 4: Dominant Pluto lines in h 
their common logarithms. 



k and p - q. Amplitudes are represented by 



Line 


Identity 


Period 

yr 


Frequency 
day- 1 


h-k 


P-< 


S 


2*pl-g8 


1.37e8 


-2.0e-ll 


-1.36 


-2.42 


W 


p2 


2.75e7 


-9.95e-ll 


-1.88 


-2.86 


R 


Pi 


3.69e6 


-7.413e-10 


-0.63 


-0.86 


T 


g8 


1.87e6 


-1.462e-9 


-1.47 


-2.63 


U 


2*g8-pl 


1.25e6 


?2.184e-9 


-2.41 


-3.30 


V 


3*g8-2*pl 


9.41e5 


-2.904e-9 


-3.05 


-3.57 



Note: The ? in the frequency of line U indicates that the sign of the frequency could 
not be determined because the phases were ill defined. 



Table 5: Pluto frequencies in e, t, and oj. Amplitudes are represented by their common 



logarithms. 












Line Identity 


Period 

yr 


Frequency 
day -1 


e 


t 


U) 


L p2-2*pl+g8 


3.44e7 


7.94e-ll 


-3.4 


-2.7 


-2.2 


M pl-p2 


4.27e6 


6.418e-10 


-2.3 


-2.4 


-1.0 


N g8-pl 


3 . 796e6 


7.2116e-10 


-1.7 


-1.8 


-.42 


3*pl-p2-2*g8 


3.419e6 


8.008e-10 


-2.9 


-3.0 


-1.7 
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Appendix 

In this appendix we supply invariable-frame center-of-mass coordinates and velocities 
of the outer planets listed at approximately 11 Myr intervals (4.0 x 10 9 days) for 
approximately 110 Myr. For each epoch the x, y, and z coordinates of the position 
and the velocity of each body are given. Positions are in AU, velocities are in AU/day. 
Caution: This table is printed with a large number of significant places. This is not 
intended to indicate the accuracy of these numbers. After 110 Myr the integration 
error in the position of Jupiter is a large fraction of an orbit. Many places are offered 
solely for the purpose of precise reproduction and testing of our results. 
We integrate Pluto as a zero-mass test particle, for the other planets we use the CHO 
masses: 



Body 


Reciprocal Mass 


Sun 


1.0 


Jupiter 


1047.355 


Saturn 


3501.6 


Uranus 


22869. 


Neptune 


19314. 


Pluto 


zero mass 



The sum of the masses of the inner planets is 0.597682 x 10~ 5 . 

The Gaussian constant is 0.01720209895. 

Our initial conditions are for J.D. 2430000.5, copied from CHO: 

J.D. 2430000.5 
Sun 

-4 06428667034226E-003 -6 . 08813756435987E-003 - 1 . 66162304226834E-006 

6 . 69048890636161E-006 -6 . 33922479683693E-006 -3 . 13202146B90767E-009 
Jnpiter 

3 . 40546614227466E+000 3 . 62978 190075864E+000 3 . 42386261766677E-002 

-5 . 697979693 10664E- 003 6 . 6181B399480116E-003 -2 . 66711392865691E-006 
Saturn 

6 . 60801 B64403466E+000 6 . 38084674686064E+000 - 1 . 36146963724642E-001 

-4.17354020307064E-003 3. 997237617481 16E-003 1 .67206320671441E-006 
Uranus 

1 . 11636331406697E+001 1 . 60373479067266E+001 3 . 61783279369968E-001 

-3.268848061S1064E-003 2.06438412906916E-003 -2. 17699042180569E-006 
Neptune 

-3.01777243405203E+001 1.91166314998064E+000 -1 .63887696621042E-001 

-2.17471786046638E-004 -3.11361111026884E-003 3.58344706491441E-006 
Pinto 

-2 . 13868977631673E+001 3 . 20719104739886E+001 2 . 49246689656096E+000 

-1 .76936677262484E-003 -2.06720938381724E-003 6.68091931493844E-004 
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Sun 

3.48223323133294E-O03 
6 . 99036269676260E-006 
Jupiter 

1 .61746471832066E+000 
-7.08766744678321E-003 
Saturn 

-1 .00147212070866E+001 
-1 .06660965907352E-003 
Uranus 

-1 .40377336981219E+001 
2 . 48032 16428B927E-003 
Neptune 

-2 .99890268409148E+001 
-6.22984243478207E-004 
Pluto 

4.46173202040767E+001 
1 .37872077689202E-004 



J.D. 4002430000.6 

•6 . 14286882438496E-003 
•6.68635026091231E-007 



6 . 42296997953432E-005 
• 1 . 84571827079439E-008 



4.96217444761191E+000 -1.88332403666667E-002 

2 . B20106B8679173E-003 3 . 534630B2614887E-006 

2 . 3482998289B446E+000 - 1 . 04707804806302E-001 

-6 . 03824327277739E-003 -6 . 436939472B8976E-006 



-1 .21373363763877E+001 
•3.21481640139462E-003 



3 . 13477343821 166E-002 
7 . 6B66370327BB8BE-O0B 



6 . 306363413204B4E+000 -3 . 4217B281320721E-001 
■3.0637000537262BE-003 -4.93621278691276E-006 

4.64963741941286E+000 -1 .11244716244026E+001 
2 . 277B6980086148E-003 -3 . 7178B740399728E-004 



Sun 

1 . 66869198406B22E-004 
B . 39160393766469E-006 
Jupiter 

-S.44729639898946E-001 
-7 . 16676004968646E-003 
Saturn 

3 . 64B64162B80684E+000 
B.32043382196174E-003 
Uranus 

1 .96630741804694E+001 
-6 . 60823674708960E-004 
Neptune 

-2 . 90434891624838E+001 
- 7 . 706008 1 4382826E-O04 
Pluto 

1 .761B2664330301E+O01 
2 .96906833418316E-003 



J.D. 8002430000. B 

-3 . 18403738476778E-003 
3.43231236267224E-008 

E . 41799893231963E+000 
-B.94418786410186E-004 

-8 . 74621730690643E+000 
1 .84433478081 138E-003 

2 . 71780478978726E+000 
3 . 760132276SB692E-003 

7 . 63163376490776E+000 
-3 . 04993242063366E-003 

-2.826600621069B1E+001 
9 . 61086077046478E-004 



1 .08836667678877E-006 
7 . 3BSS6929980973E-008 

1.13960632794B64E-002 
-6.0804787962B361E-006 

2.0613B684103364E-002 
-8 . 192677B9180138E-006 

- 2 . 93842766872691E-001 
4 . 0279904B696946E-006 

-2.86883308666771E-001 
-3 . 63038934787746E-00B 

2.07127824781669E+000 
-7 . 391616496247B0E-004 
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Sun 

2 . 16760693763746E-003 
7.61620042967966E-006 

Jupiter 

-1.91624116724281E+000 
-6.75468231460021E-003 

Saturn 

6.76826626624820E+000 
-3 . 82866204482741E-003 

Uranus 

-1.94371685964776E+001 
-1 . 44364395609248E-004 

Neptune 

-2 . 74440706887983E+001 
-1 .27960784828791E-003 

Pluto 

-1 . 20764756309941E+001 
-2 . 20618B57580004E-003 



J.D. 12002430000.6 

-7 . 22222306693326E-003 
2 . 10694058046252E-006 

4.90690123135316E+000 
-3 . 10003062639296E-003 

6 . 62826051289616E+000 
4.10289639223501E-003 

7.12459915768657E-001 
■3 . 87282609634166E-003 

1 . 18415484883449E+001 
-2 . 88666664819439E-003 

3 . 86602941387661E+001 
1 . 27298474290349E-003 



-3.1 6449640367761E-006 
-2 .63984886742416E-008 

-1.21294100709737E-002 
4 . 40061982227406E-006 

1 .36821733133484E-001 
-3 . 84208074824303E-005 

-6.67184312770713E-002 
-7.26422112586187E-006 

1 .26168495620B80E-001 
-2 . 83733014319321E-006 

-7 . 61166871262173E+000 
6 . 96668670088588E-004 



Sun 

6 . 50963831 183654E-003 
5 . 63227733444709E-006 
Jupiter 

-3 . 74866453594B65E+000 
-B . 13292664282663E-003 
Saturn 

-8.21572601979340E+000 
-2 . 14461713364704E-003 
Uranus 

1 69016814347972E+001 
1.6B074386127937E-003 
Neptune 

-2 . B6673979766374E+001 
-1 .676B6082030944E-003 
Pluto 

4 . 40762063188756E+001 
2 . 16842046624218E-004 



J.D. 16002430000.6 

-4 . 78476579136663E-003 
7 . 13351656818888E-006 



-4 . 19212523666333E-006 
3 . 47104548933087E-008 



3.34109678830964E+000 -1 .93001386469316E-002 
•5.86304160966797E-003 -B.20870372647430E-006 



3 . 67351862497295E+000 
-6.48880502566166E-003 

-6 . 56637329418294E+000 
3 . 86196567636484E-003 



9.12765874981280E-002 
6 . 49049636096309E-005 

3 . 20323467366284E-001 
-1 . 16686367761641E-006 



1 . 60830378404650E+001 3 . 91 697006367944E-001 
2.63696428483002E-003 -2.87482304637816E-006 

4.24164281996522E+000 -6.14997629409786E+000 
2 . 30927079918678E-003 -6 . 46238501664341E-004 
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Sun 

5 . 52370658890758E-003 
1 .60116973763634E-006 
Jupiter 

-4 . 68031337944689E+000 
-2.94440784864176E-003 
Saturn 

2 . 60191641873676E+000 
6. 02899732471 168E-003 
Uranus 

-1 . 19983266302817E+001 
-2 .97637197E82829E-003 
Neptune 

- 2 . 40433361 247367E+001 
-1 .86403697610088E-003 
Pinto 

5 . 11389786669684E+000 
3. 387138297981 16E-003 



J.D. 20002430000.6 

-6.86828041842782E-004 
6.69992634492122E-006 

1 .827264271 88393E+000 
-7.21603288617162E-003 

-9 . 81704969722601E+000 
1 .46797094688450E-003 

1 .66889069176537E+001 
-2 . 16413786466044E-003 

1 .77539929886424E+001 
•2.56670198246339E-003 

■2 .96862108067399E+001 
9 . 15467760686447E-006 



-1 .70676698693163E-006 
-4.13311417283429E-008 

1 .42744898234868E-002 
9.60218682533074E-006 

-8 . 68646694665737E-003 
9.84269561882108E-006 

-3 . 28798096836417E-002 
6 . 70903408470076E-006 

1.416441E0864033E-001 
2 . 16476730226097E-006 

8 . 28422936812740E+000 
-1 .37072741219186E-004 



Son 

4.23664486163340E-003 
1 .04376619239960E-006 
Jupiter 

-6.32661369019271E+000 
1 .91013028503073E-006 
Saturn 

6 . 57076550E76622E+000 
-3 . 96466139790090E-003 
Uranus 

3 . 42362483017968E+000 
3 . 9671341 1447061E-003 
Neptune 

-2 . 27646827694062E+001 
-2.04910439860953E-003 
Pluto 

-4 . 03030206035373E+000 
-2.189E2962349382E-003 



J.D. 24002430000.5 

1.88623421298253E-003 9.64322707746441E-006 
6.98832048345674E-006 -3.68770969502140E-009 



-2.93672213974529E-001 
-7 . 34692249237661E-003 



-4 . 29868168175026E-002 
-1.20668232770601E-006 



6 . 9267249606861 9E+000 - 1 . 031 67263666313E-001 

3.91677086101041E-003 4.78674667686964E-006 

•1.83666500192271E+001 -3.16498468287528E-001 

7.46490274236667E-004 -3.30616613730046E-006 



1 . 9137374644791 2E+001 
■2 .42379178407079E-003 



-2.33626787235126E-001 
3 . 23273419266869E-006 



4 . 56782699971212E+001 -1 . 39212971213227E+001 
■2.46591037673138E-004 1 .68666647969970E-006 
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San 

8.61197183605486E-003 
-1 .00069705 146846E-006 
Jupiter 

-6.27028919471399E+000 
1 .89769344379601E-003 
Saturn 

-9 . 19626163474942E+000 
-1 .78889464915172E-003 
Uranus 

8.69431097112419E-001 
-3 . 99280596022371E-003 
Neptune 

-1 .91637796884696E+001 
-2.42620431301626E-003 
Pluto 

2.37862118021672E+001 
-1 .731234B2178467E-003 



J.D. 28002430000. B 

-1.62191606336924E-003 2. 12686089274339E-006 

8 . 23053690598849E-006 3 . 22486369692167E-008 

- 1 . 3171 2395740664E+000 3 . 60776423487000E-002 

-6. 946684239091 19E-003 -1.71911380994973E-006 

2. 6036997065 1620E+000 -1 32824607238309E-001 

-5 . 26201677273889E-003 -4 . 39064061916997E-006 

1 .88607814395373E+001 -6.80667744546339E-003 

1.22497797829261E-004 -7.61186636206692E-006 

2. 3397167 4397009E+001 -3.20117401627937E-001 

-1.96325216017069E-003 -2.17896266864682E-007 

3 . 33316621708726E+001 -7 . 67648286684333E+000 

1 . 80924976864970E-003 -6 . 67807229934771E-004 



Sun 

5.29517287415151E-003 
-6 43467628626935E-006 
Jupiter 

-4.63094682765673E+000 
3 . 89595304220610E-003 
Saturn 

3.42671456683762E-001 
6 . 94415364763047E-003 
Uranus 

-6 . 93630690949669E+000 
3 . 63469430286160E-003 
Neptune 

-1 .65936849988451E+001 
-2 . 65201 199183682E-003 
Pluto 

-3 . 27389770115819E+001 
1 .95616788080668E-003 



J.D. 32002430000.6 

4.60016974301791E-003 2.77183378213466E-005 

6.11118866899664E-006 1 .09626231979361E-008 

-2 . 6961 1647359337E+000 -3 . 40764660633709E-002 

-6 . 40679989069352E-003 1 . 76809675487063E-006 

-8 . 88440201648099E+000 - 1 . 16080876733686E-002 

6.07707314116164E-004 -9.02260227034660E-006 

-1 .91426604330646E+001 3.31078321460774E-001 

■1 .27662427623986E-003 2.58817114609627E-006 

2.61231437761560E+001 -1 .22560187969806E-001 

•1.60914612769619E-003 -4.23093606662130E-006 

-2.08180336667930E+001 1 .03183956468066E+000 

• 1 . 82862231 122927E-003 7 . 44436957672349E-004 
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Sun 

1.100781401B1875E-003 
-4 . 32406444740262E-006 
Jupiter 

-3.57482684181278E+000 
5 .66517106961 183E-003 
Saturn 

8.47989O3166B1O1E+0O0 
-2 .91879138364036E-003 
Uranus 

1.51072149166B28E+001 
-2 . 62004191424286E-003 
Neptune 

-1 .48699767282820E+001 
-2 .72679B00617084E-003 
Pluto 

- 3 . 86424876637777E+001 
-1 .36226740673616E-003 



J.D. 36002430000.6 

3 . 21060447564082E-006 
3 . 81882603676273E-006 



•6.11 493629868046E- 006 
3.76961684870714E-008 



-3 . 61331816946287E+000 2 . 06248440092848E-002 
-B . 389774B6078203E-003 - 2 . 33897660976661E-006 



B.26572264021647E+000 
4 .44601096405309E-003 

1 .07814360862296E+001 
3 . 18266868462124E-003 

2.60124717866682E+001 
•1 .67686341901427E-003 



1 .34774368626477E-001 
-6 . 766771 27079669E-006 

•5.99164409177869E-002 
7 . 21662908916507E-006 

1 .097B6641971781E-001 
3 . 76693469208698E-006 



2.81416712262606E+001 -1 .26276918036680E+001 
1 . 6248241 180B361E-003 1 . 19699042013392E-004 



Sun 

6 . 1 1 103294604066E-003 
-6 -92783401638197E-006 
Jupiter 

-1 .5997B62B403686E+000 
7 . 24902527B95061E-003 
Saturn 

-7 . 80989146903828E+000 
-3.229626021323B2E-003 
Uranus 

-1 .70680418064660E+001 
1 .81263346708702E-003 
Neptune 

-1.1731181B271961E+001 
-2 .90342966186623E-003 
Pluto 

2 .99166240166969E+001 
6.92060814781919E-006 



J.D. 40002430000.6 

1.87723634619612E-003 
4 . 18682312303781E-006 

-4.76966161643362E+000 
-2 .84299406489444E-003 

6 .77368422097477E+000 
-4 . 41342438076074E-003 

■9 . 27396066298768E+000 
-3 . 43384641827910E-003 

2 .76773207866022E+001 
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